﻿using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Diagnostics;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using System.Windows.Forms;
using TaskBuilder;

namespace UniverseSim
{
    public partial class UniverseSim : Form
    {
        private const int W = 500;
        private const int H = 500;
        private const double G = 1.0E-13;
        private const double MX = 1.0;
        private const double MY = 1.0;
        private const double MDX2 = 1.0E-6;
        private const double MDY2 = 1.0E-6;

        private struct Vector
        {
            public double X;
            public double Y;

            public Vector(double x, double y)
            {
                X = x;
                Y = y;
            }
        }

        private struct Particle
        {
            public Vector A;
            public Vector V;
            public Vector P;

            public Particle(double x, double y)
            {
                A = new Vector(0.0, 0.0);
                V = new Vector(0.0, 0.0);
                P = new Vector(x, y);
            }

            public Particle(double px, double py, double vx, double vy, double ax, double ay)
            {
                A = new Vector(ax, ay);
                V = new Vector(vx, vy);
                P = new Vector(px, py);
            }
        }

        private class Cell
        {
            private static Random _r = new Random(100);

            public List<Particle> Particles = new List<Particle>();
            public Vector P;
            public Vector A;

            public Cell(int numParticles)
            {
                for (int i = 0; i < numParticles; i++)
                    Particles.Add(new Particle(_r.NextDouble() / 10 + 0.45, _r.NextDouble() / 10 + 0.5));
            }
        }

        private class Universe
        {
            public Cell[,] Cells;

            public Universe(int numParticles, int numCells)
            {
                int numHCells = (int)Math.Sqrt(numCells);
                int numVCells = numHCells;

                Cells = new Cell[numHCells, numVCells];

                for (int v = 0; v < numVCells; v++)
                    for (int h = 0; h < numHCells; h++)
                        Cells[h, v] = new Cell(numParticles / (numHCells * numVCells));
            }

            public IEnumerable<Cell[]> GetCells()
            {
                var cells = new Cell[Cells.Length];
                int i = 0;
                foreach (Cell cell in Cells)
                    cells[i++] = cell;

                yield return cells;
            }
        }

        public UniverseSim()
        {
            InitializeComponent();
        }

        private Cell[] SelectCells(Tuple<Cell[], Cell[]> cells)
        {
            return cells.Item1 ?? cells.Item2;
        }

        private Cell ProcessLocal(Cell cell)
        {
            Cell newCell = new Cell(0);

            double spx = 0.0;
            double spy = 0.0;

            List<Particle> ps = cell.Particles;
            for (int i1 = 0; i1 < ps.Count; i1++)
            {
                spx += ps[i1].P.X;
                spy += ps[i1].P.Y;

                double ax = 0.0;
                double ay = 0.0;
                for (int i2 = 0; i2 < ps.Count; i2++)
                {
                    // distance between points
                    double dx = ps[i2].P.X - ps[i1].P.X;
                    double dx2 = dx * dx;
                    if (dx2 < MDX2) dx2 = MDX2;
                    double dy = ps[i2].P.Y - ps[i1].P.Y;
                    double dy2 = dy * dy;
                    if (dy2 < MDY2) dy2 = MDY2;

                    // acceleration
                    ax += G / dx2 * Math.Sign(dx);
                    Debug.Assert(!double.IsNaN(ax));
                    ay += G / dy2 * Math.Sign(dy);
                    Debug.Assert(!double.IsNaN(ay));
                }

                // velocity
                double vx = ps[i1].V.X + ax;
                Debug.Assert(!double.IsNaN(vx));
                double vy = ps[i1].V.Y + ay;
                Debug.Assert(!double.IsNaN(vy));

                // position
                double px = ps[i1].P.X + vx;
                Debug.Assert(!double.IsNaN(px));
                double py = ps[i1].P.Y + vy;
                Debug.Assert(!double.IsNaN(py));

                // wrap
                if (px > MX || px < 0) px -= Math.Floor(px / MX) * MX;
                if (py > MY || py < 0) py -= Math.Floor(py / MY) * MY;

                newCell.Particles.Add(new Particle(px, py, vx, vy, ax, ay));
            }

            if (newCell.Particles.Count > 0)
                newCell.P = new Vector(spx / newCell.Particles.Count, spy / newCell.Particles.Count);
            return newCell;
        }

        private Cell[] ProcessGlobal(Cell[] cells)
        {
            Cell[] newCells = new Cell[cells.Length];

            for (int i1 = 0; i1 < cells.Length; i1++)
            {
                Cell newCell = new Cell(0);
                double ax = 0.0;
                double ay = 0.0;
                for (int i2 = 0; i2 < cells.Length; i2++)
                {
                    if (i1 == i2)
                        continue;

                    // distance between points
                    double dx = cells[i2].P.X - cells[i1].P.X;
                    double dx2 = dx * dx;
                    if (dx2 < MDX2) dx2 = MDX2;
                    double dy = cells[i2].P.Y - cells[i1].P.Y;
                    double dy2 = dy * dy;
                    if (dy2 < MDY2) dy2 = MDY2;

                    // acceleration
                    ax += G * cells[i1].Particles.Count * cells[i2].Particles.Count / dx2 * Math.Sign(dx);
                    Debug.Assert(!double.IsNaN(ax));
                    ay += G * cells[i1].Particles.Count * cells[i2].Particles.Count / dy2 * Math.Sign(dy);
                    Debug.Assert(!double.IsNaN(ay));
                }

                foreach (Particle p in cells[i1].Particles)
                    newCell.Particles.Add(new Particle(
                        p.P.X, p.P.Y, p.V.X, p.V.Y, p.A.X + ax, p.A.Y + ay));
                newCell.A = new Vector(cells[i1].A.X, cells[i1].A.Y);
                newCell.P = new Vector(cells[i1].P.X, cells[i1].P.Y);
                newCells[i1] = newCell;
            }

            return newCells;
        }

        private Cell[] Reallocate(Cell[] cells)
        {
            Cell[] newCells = new Cell[cells.Length];
            for (int i = 0; i < newCells.Length; i++)
            {
                newCells[i] = new Cell(0);
                int rc = (int)Math.Sqrt(cells.Length);
                int r = i / rc;
                int c = i % rc;
                double rcp1 = rc + 1.0;
                double minX = c == 0 ? double.MinValue : (1.0 * c) / rcp1;
                double maxX = c == rc - 1 ? double.MaxValue : (1.0 * c + 1.0) / rcp1;
                double minY = r == 0 ? double.MinValue : (1.0 * r) / rcp1;
                double maxY = r == rc - 1 ? double.MaxValue : (1.0 * r + 1.0) / rcp1;

                foreach (Cell cell in cells)
                    foreach (Particle p in cell.Particles)
                        if (p.P.X >= minX && p.P.X < maxX && p.P.Y > minY && p.P.Y < maxY)
                            newCells[i].Particles.Add(new Particle(
                                p.P.X, p.P.Y, p.V.X, p.V.Y, p.A.X, p.A.Y));
            }
            return newCells;
        }

        Bitmap _img = new Bitmap(W, H);
        int _i = 0;
        private void Pixellate(Cell[] cells)
        {
            if (InvokeRequired)
            {
                if (Interlocked.Increment(ref _i) % 10 == 1)
                    Invoke(new Action(() => Pixellate(cells)));
                return;
            }

            _img = new Bitmap(W, H);
            Color col = Color.Black;
            foreach (Cell c in cells)
                foreach (Particle p in c.Particles)
                    _img.SetPixel((int)(W * p.P.X), (int)(H * p.P.Y), col);

            Refresh();
        }

        private void UniverseSim_Load(object sender, EventArgs e)
        {
            Width = W;
            Height = H;

            var universe = new Universe(1000, 100);

            using (TB.Scope())
            {
                var realloc = TB.Point<Cell[]>();
                var cells = TB.Any(TB.Func(universe.GetCells()), realloc);
                var scells = TB.Branch(TB.Func(SelectCells, cells));
                var local = TB.Func(ProcessLocal, universe.Cells.Length, scells);
                var global = TB.Func(ProcessGlobal, local);
                realloc.Attach(TB.Func(Reallocate, global));

                TB.Dispatch(Pixellate, scells);
            }
        }

        private void UniverseSim_Paint(object sender, PaintEventArgs e)
        {
            Graphics g = e.Graphics;
            g.DrawImage(_img, new Point(0, 0));
        }
    }
}
